
# Template code for plot of time series data with distribution of imputated values.
# jH 16/10/07


stub<-"h"          # This is the stub name to your imputed datasets
nimp<-100          # This is the number of imputed datasets
alpha<-0.95        # This partially controls the confidence intervals.  Note .95 actually creates 90 percent CI's (Takes .05 off each side of the distribution).  
                   # If changed, you have to set this in the ci.upper and ci.lower functions also.
jitter<-0.2        # This is a graphics parameter


impdata<-read.csv("baummerge2.csv")   # This is the raw unimputed dataset.




#subcty<-unique(impdata$countrycode)

subcty<-c("MEX","CIV","TZA","ARE")

countrynames<-c("Country A","Country B","Country C","Country D","Country E","Country F")


# These two functions are no longer used

ci.upper<-function(data){
  alpha<-0.95   #local alpha, not global alpha, makes apply easier, and yet cumbersome
  a<-sort(data)
  b<-a[round(length(a)*alpha)] 
  return(b)
}

ci.lower<-function(data){
  alpha<-0.95
  a<-sort(data)
  b<-a[round(length(a)*(1-alpha))] 
  return(b)
}



period<-c(1959,2000)


yearseq<-(period[1]+1):(period[2]-1)
xyears<-rbind(yearseq,yearseq,NA)
delta<-0.1
xyearsdelta<-rbind(yearseq-delta,yearseq+delta,NA)

l<-1
w<-4
lw<-l*w
par(mfrow=c(l,w))   # Creates a window of l-by-w subwindows


for(i in 1:length(subcty)){

  #allmain<-countrynames[i]

  flag<-impdata$countrycode==subcty[i] & impdata$year>period[1] & impdata$year<period[2]
  tempdata<-impdata[flag, ]
  allmain<-tempdata[1,"countryname"]

  lowerb<-min(tempdata$femlifeex,na.rm=TRUE) - 0.1*abs(min(tempdata$femlifeex,na.rm=TRUE))
  upperb<-max(tempdata$femlifeex,na.rm=TRUE)*1.1

  allylim<-c(lowerb,upperb)
 
  r1<-is.na(tempdata$femlifeex)
  impeddata<-matrix(0,period[2]-period[1]-1,nimp)
  stackdata<-matrix(0,period[2]-period[1]-1,nimp)


  for(j in 1:nimp){
    temp.imp<-read.csv(file=paste(stub,j,".csv",sep=""))
    impflag<-temp.imp$countrycode==subcty[i] & temp.imp$year>period[1] & temp.imp$year<period[2]
    impeddata[,j]<-temp.imp$femlifeex[impflag]
  }

  for(i in 1:nrow(stackdata)){
    stackdata[i,]<-sort(impeddata[i,])
  }

#  s.upper <-apply(impeddata,1,ci.upper)
#  s.median<-apply(impeddata,1,median)
#  s.lower <-apply(impeddata,1,ci.lower)

   cilr<-c(round(nimp*(1-alpha)),round(nimp*alpha))

   missing<-is.na(tempdata$femlifeex)
   years.missing<-xyears[1,missing]

   s.upper <-  stackdata[,cilr[1]]
   s.median <- stackdata[,round(nimp/2)]
   s.lower <-  stackdata[,cilr[2]]

# If you uncomment the piece below, it shows all the imputed values.

#  # Imputed Values   
#  plot((years.missing*matrix(1,length(years.missing),cilr[1]-1))+runif((length(years.missing))*(cilr[1]-1),-jitter,jitter),stackdata[missing,1:(cilr[1]-1)],main=allmain,xlab="Year",ylab="Female Life Expectancy",ylim=allylim,xlim=period,type="p",col="salmon")
#  par(new=TRUE)
#  plot((years.missing*matrix(1,length(years.missing),nimp-cilr[2]))+runif((length(years.missing))*(nimp-cilr[2]),-jitter,jitter),stackdata[missing,(cilr[2]+1):nimp],main=allmain,xlab="Year",ylab="Female Life Expectancy",ylim=allylim,xlim=period,type="p",col="salmon")
#  par(new=TRUE)
#  plot((years.missing*matrix(1,length(years.missing),cilr[2]-cilr[1]+1))+runif((length(years.missing))*(cilr[2]-cilr[1]+1),-jitter,jitter),stackdata[missing,cilr[1]:cilr[2]],main=allmain,xlab="Year",ylab="Female Life Expectancy",ylim=allylim,xlim=period,type="p",col="grey")
#  par(new=TRUE)

  # Confidence Intervals  
  plot(xyears[,r1],rbind(s.upper[r1],s.lower[r1],NA),main=allmain,xlab="Year",ylab="Female Life Expectancy",ylim=allylim,xlim=period,type="l",col="slategrey")
  par(new=TRUE)
  plot(xyearsdelta[,r1],rbind(s.median[r1],s.median[r1],NA),main=allmain,xlab="Year",ylab="Female Life Expectancy",ylim=allylim,xlim=period,type="l",col="slategrey")
  par(new=TRUE)

  # Some (m) Imputed Values
  m<-5
  plot((years.missing*matrix(1,length(years.missing),m)),impeddata[missing,6:(5+m)],main=allmain,xlab="Year",ylab="Female Life Expectancy",ylim=allylim,xlim=period,type="p",col="blue")
  par(new=TRUE)

  # Observed Data
  plot(tempdata$year,tempdata$femlifeex,main=allmain,xlab="Year",ylab="Female Life Expectancy",ylim=allylim,xlim=period,pch=16)
  
# or pch=16,col="red2", used with "skyblue2" also look reasonable.


# If you've got multiple pages of graphs to look through, you can uncomment the piece below:

#  if(round(i/lw)*lw==i){
#    par(ask=TRUE)
#  }else{
#    par(ask=FALSE)
#  }

}

  
